Chemical accuracy for the van der Waals density functional 
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The non-local van der Waals density functional (vdW-DF) of Dion et al. [Phys. Rev. Lett. 
92, 246401 (2004)] is a very promising scheme for the efficient treatment of dispersion bonded 
systems. We show here that the accuracy of vdW-DF can be dramatically improved both for 
dispersion and hydrogen bonded complexes through the judicious selection of its underlying exchange 
functional. New and published exchange functionals are identified that deliver much better than 
chemical accuracy from vdW-DF for the S22 benchmark set of weakly interacting dimers and for 
water clusters. Improved performance for the adsorption of water on salt is also obtained. 



London dispersion interactions are ubiquitous in na- 
ture contributing to the binding of biomolecules such as 
DNA, molecular crystals, and molecules on surfaces. The 
accurate description of dispersion, which often occurs in 
conjunction with hydrogen bonds, is a major challenge 
for many electronic structure theories. Density functional 
theory (DFT), the most widely used electronic struc- 
ture theory, often doesn't meet this challenge. Indeed, it 
is well-established that popular generalized gradient ap- 
proximation (GGA) or hybrid exchange-correlation func- 
tionals are inadequate for the description of dispersion 
interactions. Many schemes have been developed that al- 
low dispersion to be accounted for within DFT in a more 
or less approximate manner (see, e.g. [E S S S S 0] ) • 
One of the most promising and rigorous methods is the 
non-local van der Waals density functional (vdW-DF) of 
Langreth and Lundqvist and co-workers [l[. 

In vdW-DF the non-local correlation is calculated so 
that the exchange-correlation energy takes the form 
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where E^ GA is the GGA exchange energy. In the original 
vdW-DF this is obtained with the revised version of the 
Perdew, Burke, and Ernzerhof (PBE) Q functional from 
Zhang and Yang (revPBE) |{|. E^° A accounts for the 
local correlation energy obtained within the local density 
approximation (LDA), and E™ 1 is the non-local correla- 
tion energy. The formula for E^ 1 is based on electron 
densities interacting via a model response function, the 
particular form of which is still a subject of research [l(| • 
The vdW-DF has been applied to a wide variety of sys- 
tems where dispersion is important (see ref. [11] for a 
review) and recent algorithmic developments [12j have 
made it only marginally more computationally expen- 
sive than a regular GGA. However, in many important 
circumstances the current vdW-DF is simply not accu- 
rate enough. For example, for the S22 dataset [3] (a 
set of 22 weakly interacting dimers mostly of biological 
importance) it yields a mean absolute deviation (MAD) 



of ^60 meV [14[ compared to coupled cluster reference 
data. This is outside the so-called "chemical accuracy" 
of 1 kcal/mol or ~43 meV and inferior performance to 
other DFT-based dispersion correction schemes 0, H, 0| ■ 
Water clusters, important for atmospheric chemistry and 
liquid water, are another example where vdW-DF sub- 
stantially underbinds (by ~20% compared to accurate 
reference data) and in terms of absolute dissociation en- 
ergies is worse than a regular GGA such as PBE 15. 

Recognizing that the interaction energies obtained 
with vdW-DF depend on the exchange functional incor- 
porated within it [H, Ei , 12 1 , we aimed to improve vdW- 
DF by exploring and developing alternative exchange 
functionals to the original revPBE. We take a pragmatic 
approach, we use E^ 1 in its regular form and search for 
an exchange functional that combines with it to give pre- 
cise energies for a wide range of systems. To this end 
we first use the S22 dataset since it includes a variety 
of weakly bonded dimers for which accurate interaction 
energies and structures have been established and 
so provides a tough test for molecular simulation meth- 
ods. Following this we test our methods on two complex 
systems where dispersion interactions are crucial: wa- 
ter hexamers and water adsorbed on NaCl(OOl). From 
these studies we propose three new exchange functionals, 
that when incorporated within vdW-DF offer vastly im- 
proved interaction energies compared to those from the 
original vdW-DF. The new functionals, which are easy 
to implement and come at no extra cost, make vdW-DF 
competitive with all other DFT-based methods for the 
treatment of weak interactions. We hope that this study 
lays the foundations for further improvements of vdW- 
DF and enables more accurate treatments of dispersion 
and hydrogen bonded systems, for example, liquid water 
and ice. 

Throughout, we calculate the vdW-DF energies non- 
self-consistently in two steps. First, VASP 5.2 flp , Hij 
calculations with a given exchange functional [20( and 
PBE correlation functional are performed. Second, the 
VASP electron density is used to determine the vdW 
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correction using JuNoLo 2l|. We find that the mag- 
nitude of the vdW correction is rather insensitive to the 
underlying density used [22| . Therefore, density from 
B86 (exchange) and PBE (correlation) calculations was 
used for all functionals except PBE and revPBE, where 
density from the respective exchange-correlation func- 
tional was used. Care was taken with the VASP calcu- 
lations to ensure that converged energies were obtained, 
which involved the use of hard projector-augmented wave 
(PAW) [H 0] potentials, an 800 to 1000 eV cut-off, 
dipole corrections, and 20 to 25 A 3 unit cells. Since the 
efficient self-consistent calculation of vdW-DF energies 
has only very recently become possible [IH , we checked at 
the latter stages of this study how the non-self-consistent 
and self-consistent interaction energies differ with the 
grid-based GPAW code 0. For the S22 dataset the non- 
sclf-consistent and self-consistent interaction energies are 
within 1.5 meV, except for the large dispersion bonded 
dimers (dimers 11-15) where the differences are <4 meV. 
For the water hexamers the non-self-consistent and self- 
consistent interaction energies are within ^2 meV. 

Let us first examine the results for the standard form of 
the vdW-DF where revPBE exchange is used. Through- 
out, we denote a combination of an exchange functional X 
with vdW correlation as X-vdW, hence we refer to vdW- 
DF as revPBE-vdW. The differences in the revPBE-vdW 
and reference interaction energies for each of the dimers 
in the S22 dataset are shown in Fig. [TJ One can see from 
Fig. [1] that with revPBE-vdW most of the dimers are 
substantially underbound. The MAD is 65 meV (Ta- 
ble H|, which is in good agreement with the MAD of 
60 meV in Gulans et al. [14J. The errors for the indi- 
vidual hydrogen bonded (HB), dispersion bonded (DB), 
and mixed dispersion and hydrogen bonded (MB) sub- 
sets are all quite large at 106, 52, and 38 meV, re- 
spectively. Further, this functional yields a very large 
"Range" of errors (i.e., the difference between the largest 
and smallest errors) of 163 meV and does not provide a 
good balance between H bonding and dispersion. This is 
shown by A H ~ D which gives the difference between the 
mean deviations (MD) of the HB and DB subsets (i.e., 
A H ~ D =MD(HB)— MD(DB)). The smaller this quantity, 
the better the balance between the different types of 
bonding. The rather large value of A H_D for revPBE- 
vdW of 57 meV reveals that on average the HB dimers 
are underbound compared to the DB dimers. Therefore, 
the original vdW-DF does not deliver chemical accuracy 
for either systems held by dispersion or H bonds, it yields 
a large range of errors, and on average underbinds the H 
bonded compared to the dispersion bonded systems. 

We now consider alternatives to revPBE for E^ GA in 
([1]), discussing the results from just a few of the most in- 
teresting functionals, namely PBE, Becke86 (B86) [26j ]. 
and Becke88 (B88) As can be seen from Table Hall 
functionals yield smaller MADs than revPBE-vdW and 
overall the performance is qualitatively different. For ex- 



TABLE I: Mean absoiute deviations from the reference 
data [13j for the S22 set for vdW-DF with various exchange 
functionals ("Method"). MADs are given for the whole set 
("MAD Total"), the hydrogen bonded ("MAD HB"), dis- 
persion bonded ("MAD DB"), and mixed bonding subsets 
("MAD MB"). For each functional we also report the dif- 
ference between the largest and smallest deviations from the 
reference data ("Range") and the difference in the mean de- 
viations of hydrogen bonded and dispersion bonded subsets 
(A H ~ D ). The new functionals introduced in this study are in 
the last three rows. All values are in meV. 
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revPBE-vdW 


65 


106 


52 


38 


163 
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B88-vdW 


62 


76 


61 


48 


124 


16 


PBE-vdW 


54 


33 


83 


42 


117 


50 


B86-vdW 
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10 


112 
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PBE«.=l-vdW 


21 


35 


19 


10 


102 


54 


optPBE-vdW 


15 


21 


16 
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77 


37 


optB88-vdW 


10 


13 


10 


7 


44 
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ample, in contrast to revPBE-vdW, PBE-vdW systemat- 
ically overbinds the dimers. Although the MAD of B88- 
vdW is only marginally smaller (62 meV) than that of 
revPBE-vdW, the A H_D is reduced from 57 to 16 meV. 
And of most interest, B86-vdW, yields an overall MAD 
of just 23 meV. This is the lowest MAD obtained from 
all published functionals considered and a substantial im- 
provement over revPBE. 
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FIG. 1: Differences in interaction energies for vdW-DF 
(£7int[DFT]) with various exchange functionals from the 
CCSD(T) reference data (£i„ t [ACCSD(T)]) Qj|. We show 
data for revPBE, B88, PBE, B86 and three new exchange 
functionals "PBEk=1", "optPBE", and "optB88". 

Can the errors on the S22 set be further reduced? The 
contrasting performance of PBE-vdW and revPBE- vdW, 
which was also observed by Gulans et al. [1J|, provides 
the necessary physical insight to identify improved ex- 
change functionals. As we know, the GGA exchange en- 
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ergy density is given by e x (n, s) = e x (n)F x (s), where 
e x DA ( n ) i s the LDA exchange energy density and F x (s) 
is the enhancement factor that depends on the reduced 
density gradient ,s = | Vn|/2(37r 2 ) 1 / 3 n 4 / 3 . The enhance- 
ment factors of PBE and rcvPBE have the same form: 



2.5 r 
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The parameter /i is also the same and so the function- 
als differ only in the value of the parameter k. revPBE 
has a larger value of k than PBE (V cvPBE =1.245, 
k pbe =0.804), which causes F x to rise more rapidly with 
rcvPBE than PBE (see Fig. [5] to see the enhancement 
factors). As a consequence, regions with large reduced 
density gradients are stabilised more with revPBE than 
PBE, which in turn leads to weaker interactions with 



rcvPBE (see ref. [28[ for a more detailed discussion on 
this issue). Therefore, in principle, a simple strategy 
for obtaining improved interaction energies is to iden- 
tify an exchange functional intermediate between PBE 
and revPBE. To this end we varied k from the PBE to 
revPBE values (in 0.05 increments) and calculated in- 
teraction energies within vdW-DF for the complete S22 
dataset. A value of k=1.00 resulted in the smallest MAD 
of only 21 meV. We dub this new exchange functional 
"PBEk=1" [5g|]. Pushing the PBE-style (i.e., PBE and 
its various revised forms) functionals yet further we var- 
ied fi and k, and also considered other forms of the en- 
hancement factor. After optimization we obtained an ex- 
change functional "optPBE" that yielded a MAD of only 
15 meV. This functional turned out to be a 95% PBE 
and 5% RPBE @ combination with /i = 0.176, and 
k = 1.05. The enhancement factor is shown in Fig.O 

The two PBE-style functionals introduced above of- 
fer substantial improvements over revPBE. However, 
they still exhibit large errors in A H ~ D and overbind 
the methane dimer (dimer 8, by 25 meV or 108% with 
PBEk=1 and 19 meV or 85% with optPBE). Since B88 
is free from these deficiencies we explored optimised ver- 
sions of it. The B88 exchange enhancement factor can 
be written as 



F® 88 {s) = 1 + /is 2 /(l + /3s arcsinh(cs)) 



(3) 



where c = 2 4 / 3 (3tt 2 ) 1 / 3 , /i £3 0.2743, and = 
9/i(6/7r) 1 / 3 /(2c). As B88 underbinds the dimers, we 
modified the ratio fi/P to lead to increased binding, re- 
sulting in an optimal fj,/(3 of 1.2 and a fi of 0.22. The 
new exchange functional, which we dub "optB88" , yields 
a MAD of only 10 meV, an accurate binding energy for 
the methane dimer, and similar mean deviations for all 
three subsets. Of all the functionals considered, optB88 
is the most accurate for the S22 dataset. 

We now consider whether the improved performance 
of vdW-DF carries over to other systems and start with 
gas phase water hexamers. Hexamers are interesting for 
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FIG. 2: Enhancement factors of some of the exchange func- 
tionals discussed. PBE, revPBE, and the new PBEk=1 all 
have the form given by ([2} but differ in the value of k. optPBE 
is a combination of PBE and RPBE with parameters opti- 
mized for the S22 dataset. Similarly optB88 has the form 
given by (J5J and is again parameterized for the S22 dataset. 
The enhancement factor of B86 is very similar to PBEk—1 
and is not shown for clarity. 



many reasons (e.g. they are key constituents of the con- 
densed phases of water) and not least because they pro- 
vide a tough test for DFT [l6| . In particular many func- 
tionals (including those widely used to study water such 
as PBE and BLYP) incorrectly predict that a "cyclic" 
or "book" cluster has the lowest energy, in contrast to 
coupled cluster which favors a "prism" structure. The 
energy differences between the various isomers are very 
small (< 10 meV) and only when dispersion is taken into 
account is the correct energy ordering recovered (lr3 | . We 
have tested the various functionals discussed above on the 
water hexamers and the results are reported in Table HT1 
Compared to pure PBE the vdW-DF improves the rela- 
tive energies of the hexamers and B86-vdW and the three 
new functionals yield considerably improved absolute dis- 
sociation energies over revPBE- vdW 15]. In particular 
the dissociation energies for optPBE- vdW are essentially 
identical to those obtained with ACCSD(T). The fact 
that there is now a functional which predicts both accu- 
rate absolute and relative energies for water hexamers is 
very encouraging and makes optPBE-vdW an interest- 
ing prospect for condensed phase simulations of water. 
Finally we note that for one of the new exchange func- 
tionals (optB88-vdW) we have performed self-consistent 
geometry optimizations which resulted in only slightly 
different dissociation energies (net differences of ~7 meV, 
Table II) and, moreover, very similar geometries [30j ]. 

We have also applied the new functionals to another 
important class of problem, namely adsorption on sur- 
faces. The accurate determination of adsorption ener- 
gies is an issue of central importance to many disciplines. 
However, in general, there is a paucity of accurate refer- 
ence data. Water on NaCl(OOl) is an exception where an 
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TABLE II: Dissociation energies (meV/HbO) for four low 
energy isomers of the water hexamer calculated using 
ACCSD(T), PBE (with no vdW correction), and the vdW- 
DF with various exchange functionals. Unless indicated 
otherwise MP2 geometries from ref. pi| were used. The 
ACCSD(T) data in this table was computed as part of this 
study in the standard way, i.e., the MP2-CCSD(T) differ- 
ence at the triple zeta level (aug-cc-pvtz basis set) was added 
to the MP2 complete basis set dissociation energies. 
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revPBE-vdW a 


-280 
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-277 
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B88-vdW 
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PBE-vdW 


-380 


-378 


-372 


-358 


B86-vdW 


-328 


-327 


-324 


-314 


PBE K =l-vdW 


-326 


-325 


-322 


-313 


optPBE-vdW 


-335 


-334 


-332 


-323 


optB88-vdW 


-347 


-347 


-344 


-334 


optB88-vdW 6 


-352 


-354 


-349 


-339 



"From ref. 

'Optimized self-consistently with GPAW 



adsorption energy of -487±60 meV at the ACCSD(T) 
level was recently obtained using an embedded cluster 
approach [3lj ]. Using the geometry from ref. (3lT | and a 
slab model of the surface we computed adsorption ener- 
gies for vdW-DF with revPBE exchange, B86 exchange 
and the three new functionals. In contrast to revPBE- 
vdW, which yields an adsorption energy of —334 meV, 
the alternative choices of exchange predict adsorption en- 
ergies of -413 (PBEK=l-vdW) to -424 meV (optB88). 
Although the values presented are below the lower end of 
the error bar on the reference adsorption energy, they are 
closer to it than revPBE-vdW and a regular GGA such as 
PBE (—328 meV). This is an encouraging development 
with scope for improvement. 

In summary, we have shown that the accuracy obtained 
from vdW-DF for a range of systems can be greatly im- 
proved by making alternative choices for the GGA ex- 
change component. Based on a combination of physical 
insight and optimization three new exchange functionals 
have been proposed (PBEk=1, optPBE, and optB88). 
PBEk=1 is the simplest alternative, optPBE is an op- 
timized PBE-style functional that in addition to a low 
MAD on the S22 set yields very precise results for the 
water hexamers, and optB88 yields the overall best per- 
formance on the S22 dataset. We hope that this study 
lays the foundations for further improvements of vdW- 
DF and will enable more accurate treatments of a wide 
variety of dispersion and H bonded systems, such as liq- 
uid water. Along with the recent efficiency improvements 
[LH we suggest that vdW-DF is now a serious medium- 



term contender for high precision simulations before more 
rigorous approaches such as the random phase approxi- 
mation or quantum Monte Carlo become routine. 
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